gusucode.com > 阵列信号处理书的源码 > MATALB 程序/4.面阵中二维角度估计 Unitary -ESPRIT算法MATLAB程序/Unitary_esprit.m

    % Developed by xiaofei zhang (南京航空航天大学 电子工程系 张小飞）
% EMAIL:zhangxiaofei@nuaa.edu.cn

clear all
close all
derad = pi/180;
radeg = 180/pi;
twpi = 2*pi;
kelm = 8;               % 
dd = 0.5;               % 
d=-(kelm-1)/2*dd:dd:(kelm-1)/2*dd;     % 
iwave = 3;              % number of DOA
theta1 = [10 20 30];  
theta2 = [20 25 15];% DOA
snr = 20;              % input SNR (dB)
n = 200;                % 
A0=exp(j*twpi*d.'*(sin(theta1*derad).*cos(theta2*derad)))/sqrt(kelm);
A1=exp(j*twpi*d.'*(sin(theta1*derad).*sin(theta2*derad)))/sqrt(kelm);%%%% direction matrix
S=randn(iwave,n)
X0=[];
for im=1:kelm
      X0=[X0;A0*diag(A1(im,:))*S];
end
X=awgn(X0,snr,'measured');
L=iwave;
J1=eye(kelm-1,kelm);
J2=flipud(fliplr(J1));
Q=qq(kelm);
Y=kron(Q',Q')*X;
Q0=qq(kelm-1);
K1=real(Q0'*J2*Q);
K2=imag(Q0'*J2*Q);
I=eye(kelm);
Ku1=kron(I,K1);
Ku2=kron(I,K2);
Kv1=kron(K1,I);
Kv2=kron(K2,I);
E=[real(Y),imag(Y)];
Ey=E*E'/n;
[V,D]=eig(Ey);
EVAs =diag(D).';
[EVAs,I0] = sort(EVAs);
EVAs=fliplr(EVAs);
EVs=fliplr(V(:,I0));
Es=EVs(:,1:L);
fiu=pinv(Ku1*Es)*Ku2*Es;
fiv=pinv(Kv1*Es)*Kv2*Es;
F=fiu+j*fiv;
[VV,DD]=eig(F);
EVA = diag(DD).';
u=2*atan(real(EVA))/pi;
v=2*atan(imag(EVA))/pi;
theta10=asin(sqrt(u.^2+v.^2))*radeg
theta20=atan(v./u)*radeg